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INTRODUCTION 

One of the major tasks in dinical pharmacology is to estimate population charac- 
tenstics, such as the mean and inter-individual variability of pharmacokinetic (PK) 
and phannacodynarnic (PD) parameters in a patient population. Population anal- 
ysis techniques try to quantify this variability within the population and account 
for it in terms of patient specific covariates such as age, weight, disease state, con- 
current therapy and others. The derived population PK or PD model provides the 
necessary information required to optimize individual dose: for initial dose the 
population mean PK and PD parameters and their relationships with patient specific 
covariates suffice; for subsequent dose adjustment, the probability distribution of 
inter-mdividual random effects which capture inter-individual differences that are 
not explicable by patient specific covariates; and the residual variability due to mea- 
surement errors and intra-individual variability can be used, along with observed 
individual responses to the initial dose, to arrive at an adjusted dose. 

Population analysis according to (nonlinear) mixed effects models, as a means 
of obtaining population PK-PD parameter estimates, has gained considerable pop- 
ularity during the past decade because it permits PK-PD characterization of actual 
patient populations. It permits such characterization because it can make use of 
only a few samples from each individual collected during routine clinical care or 
as a supplement to a study designed for other purposes. Certain patient groups. 
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especially the elderly, the critically 111, and children may be hard to study under the 
rigorous conditions of well controlled clinical trials, and for them, sparse sample 
based methods may be the only practical way to determine population PK-PD. 

The discovery of a population model that adequately describes a given data 
set is often a complicated and time consuming task. This is partly due to the fact 
that the relationships of greatest interest, those between PK or PD parameters and 
covariates, are not directly observed; only the consequences of those relationships, 
concentrations or effects vs time, are observed. The modeling task consists not only 
in identifying those covariates (patient specific characteristics) that significantly in- 
fluence the PK or PD parameters but also consists in characterizing the shape of 
the relationships between covariates and parameters. In a previous paper [1] we 
presented some model-building techniques that simplify this task. In this paper we 
review that previous work, and also extend it. We present here methods that, (i) 
efficiently find and depict relationships between PK or PD parameters and patient 
covariates when they exist, (ii) screen for potential interactions between covariates, 
and (iii) provide diagnostics for detecting individuals that have an unusually large 
influence on the choice of covariates or shape of relationships. While the NONMEM 
program [2,3] is used for the work reported here, the ideas and techniques we dis- 
cuss are general and are easily transferred to other nonlinear mixed effects analysis 
platforms. 

BASIC POPULATION MODEL 

The population model has two parts: (i) a model for the response of an individual 
within the population (individual model) and (ii) a model relating the individual PK- 
PD parameters to observable patient characteristics (population model). 

The jth measurement (e.g., plasma concentration or effect) for the ith individual 
can be related to the PK or PD parameters in the following way (individual model) 

y.i = /(<k,*i;) + c O' 0) 

where / is a function (PK-PD model) describing the expected value of the response for 
a given parameter vector fa as a function of the covariate values x ijt fa is the vector of 
the model parameters (such as CI, EC50, etc.) for the ith individual and z h is the set of 
known covariates for this individual (such as dose, time, etc.). The term c,j accounts 
for the (random) error between the true value and the corresponding measurement. 
It is generally assumed that aj are independently normally distributed with either 
a constant variance a 1 for all observations (additive residual error model) or with a 
variance a 2 f (<k,x t >) 2 that is a function of the predicted drug concentration or effect 
measurement (proportional residual error model). 

An important part of variability in pharmacokinetic and pharmacodynamic ex- 
periments is due to differences between patients. The individual PK-PD parameters 
can be related to observable patient characteristics in the following way (population 
model) 

fa = g(0,z i ) + r h (2) 

where g is a known function that describes the expected value of fa as a function of 
known individual specific covariates z if such as age, weight, disease states, etc., and 
the vector of population average parameters 0. The vector r/ ( express the random 
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difference between the population predictions and individual PK-PD parameters, it 
is assumed that 7/, are independently and multivariately normally distributed with 
mean zero and variance-covariance matrix $i (diagonal elements denoted u\ where A- 
denotes the Arth parameter). Frequently, an exponential inter-individual error model 
is used which prescribes a log-normal distribution for the PK-PD parameters 



^-^^•exptr/,) (3) 
in which case a;* is an approximate coefficient of variation. 



DATA 

Data on four drugs were used to illustrate our model building and evaluation 
methods (these are the same data sets used in our previous report [1]). The analysis 
reported here for these drugs are not meant to be definitive but are presented only to 
illustrate the procedures which are the subject of this paper. 

The first drug is ibuprofen, a nonsteroidal anti-inflammatory. Ibuprofen was 
administered as a single oral dose of 5 mg/kg or 1 0 mg/kg to 92 children with febrile 
illness. Two to 6 plasma samples were obtained from each individual up to 8 hours 
post dosing (total of 411 plasma samples). The following demographic data were 
recorded: gender (SEX), male (53), female (40); race (RACE), Caucasian (17), black 
(76); location (LOC), clinic (65), hospital/inpatient (28); food since dosing (FED), no 
(25), yes (68); dose level (DRG), 5 mg/kg (45), 10 mg/kg (48); height (HT), 64-148 
cm; weight (WT), 5.6-45.5 kg; age (AGE), 3-132 months. These data were originally 
reported in [4], and were kindly supplied by Bristol Myers, Hillside, New Jersey. 

The second drug is the antiarrhythmic quinidine, administered orally to 136 
hospitalized men for various arrhythmias. Plasma samples were obtained for routine 
clinical purposes (total of 361 samples). The following demographic data were 
recorded: race (RACE), Caucasian (91), latin (35), black (10); smoking (SMO), no (91), 
yes (45); ethanol abuse (El), no ethanol abuse or social drinker (90), current ethanol 
abuse (16), history of ethanol abuse (30); congestive heart failure (HF), no or mild 
(56), moderate (40), severe (40); creatinine clearance (RF), < 50 ml/min (48), > 50 
ml/min (88); weight (WT), 41-119 kg; height (HT), 154-202 cm; age (AGE), 42-92 
years; a-l-acid glycoprotein concentration (GLP), 39-316 mg/dl. These data were 
originally reported in [51, and were kindly supplied by Thomas M. Ludden. 

The third drug is the antihypertensive prazosin, administered twice daily in a 
dose range of 1 mg per day up to 20 mg per day in patients with mild, moderate 
or severe hypertension. Each patient started with the lowest dose and the dose was 
increased until the patient's arterial pressure was controlled. After at least a period 
of eight weeks continuous therapy the PKs of prazosin were studied in 64 patients 
during a 12 hr dosing interval. After 4 weeks, 35 patients were restudied. A total 
number of 887 plasma samples were obtained. The following demographic data 
were recorded: gender (SEX), male (63), female (36); race (RACE), Caucasian (67), 
black (32); visit number (VIS), visit 15 (53), visit 16 (46); smoking (TOB), no (74), 
yes (25); prior therapy (PT), no (3), yes (96); co-treatment with hydrochlorothiazide 
(HCTZ), no (44), yes (55); co-treatment with propranolol (PROP), no (85), yes (14); 
other concomitant therapy (CON), no (24), yes (75); food (FF), fed (50), fasting (49); 
height (HT), 140-188 cm; weight (WT), 51-139 kg; age (AGE), 24-69 years; serum 
creatinine (SECR), 0.6-1.8 mg/dl. These data were- originally reported in [6], and 
were kindly supplied by Eli Lilly, Indianapolis, Indiana. 
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The last drug is the broad spectrum antibiotic pefloxacin. The PK of pefloxacin 
were studied in 74 critically ill patients with infections with organisms sensitive to 
pefloxacin. The drug was administered twice daily, 200 mg or 400 mg, in 1 hr infu- 
sions. For each individual three plasma samples were obtained during a 12 hr dosing 
interval. In total 113 dosing intervals were studied (337 plasma concentrations). The 
following demographic data were recorded: age (AGE), 18-84 years; weight (WT), 
43-125 kg; creatinine clearance (CLCR), 0.4-312 mL/min; Glasgow score (GLAS), 
3-15; simplified acute physiology score (SAPS), 1-26; albumin (ALB), 17-40 g/L; 
bilirubin (BIL), 4-150 amol/L; alanine amino transferase (ALAT), 3-200 IU/L; al- 
kaline phosphatase (AP), 32-615 IU/L; prothrombin level (PT), 36-100 % normal; 
systolic blood pressure (SPB), 60-175 mmHg; heart rate (HR), 60-170 beats/min; ar- 
tificial ventilation (AV), no (42), yes (71); gender (SEX), male (86), female (27); center 
where study was performed (CEN), center 1 (44), center 2 (69); administration of high 
dose catecholamine (CAT), no (92), yes (21). These data were originally reported in 
[7], and were kindly supplied by Rhone Poulenc Rorer, Antony, France. 

MODEL BUILDING PROCEDURE 

The population model characterizing the relationships between patient specific 
characteristics (covariates) and PK-PD parameters is derived using the stepwise 
approach previously described [1]. Briefly, in a first step individual empirical Bayes 
estimates of PK-PD parameters are obtained based on a prior NONMEM fit using no 
covariates. Subsequently, the relationship between the individual PK-PD parameter 
estimates and candidate covariates is evaluated outside NONMEM using a flexible 
semi-parametric regression model: the generalized additive model (GAM). Finally, 
the population analysis is completed using NONMEM on the basis of the GAM 
analysis of the previous step. This procedure is an extension of the graphical method 
proposed by Maitre et al. [8]. The model building procedure described above is 
demonstrated in full in this paper using ibuprofen. The results for the other drugs 
are used to address specific issues and to evaluate the robustness of the approach. 

STRUCTURAL PK-PD MODEL 

In the first step of the analysis the structural PK-PD model without any covari- 
ates is derived and fitted to the data. For ibuprofen, using NONMEM (version 4, 
level 1.0), the time profile of plasma concentrations is best described using a one- 
compartment model with first order absorption (ADVAN2, see Fig. 1). The model 
is parameterized in clearance (CI), volume of distribution (V) and the first order 
absorption rate constant ka. An exponential error model is used to characterize 
the inter-individual variability in the PK parameters. A proportional error model is 
used for the residual variability. The values of the population parameters 0, a 2 , and ft 
are estimated using the the so-called first order method in NONMEM. Subsequently, 
these parameter values are used to obtain empirical Bayes estimates of the individual 
PK parameters & using the POSTHOC step in NONMEM [9,10]. 

COVARIATE MODELS 

In the second step a regression model is derived that describes the dependence of 
the individual PK-PD parameter estimates (i.e., the elements of <fc) on the candidate 
covariates zu . . ■ zni- With the estimates of fa from the Bayes estimation step treated 
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Fig 1 Time profile of ibuprofen plasma concentrations. The solid line shows the fit of the struc- 
tural PK-model without any covariates (one-compartment mode! with first order absorption) using 
NONMEM to the actual plasma concentrations (•) obtained from 92 children with febrile illness. The 
plasma concentrations were normalized for a dose of 100 mg. 



as data, this step corresponds to the classical regression problem of variable selection 
and we take advantage of the recent work done on this problem by others [11]. 
The relationship between covariates and a specific model parameter can in the most 
general case be described by a multidimensional smooth function </(..) 

& = flr(zw...zjv;) + i?i (4) 

Without some constraints, such a general description of the data is non-identi- 
fiable, and even if identifiable, would be difficult to interpret in dimensions higher 
than two. Furthermore, the number of points required to fill a space to a given density 
grows exponentially with the dimension of the space, but usually data do not. This 
makes the estimation of complex multivariate functions using sparse noisy data 
highly inaccurate. A multiple linear regression model is an often used simplification 

N 

n=l 

However, this model makes a strong assumption, often unjustifiable on scientific 
grounds, that a linear relationship exists between parameters and covariates. A 
more general approach uses a generalized additive model (GAM) [11,12] 

N 



where g n {.) is an arbitrary univariate function. The use of splines for g n {.) is very 
appealing because splines are both flexible and convenient to use. (In this paper, 
wherever splines are mentioned, it is understood that we refer to a natural cubic 
spline with two internal knots at the 33% and 66% quantiles of the argument variable, 
equivalent to 4 degrees of freedom [1 3].) Because of the additive structure of a GAM, 
the interpretation of the analysis is straightforward, and the independent role of the 
various covariates can easily be displayed graphically. 
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Table I. Change in residual sum of squares (ARSS) when the individual (posterior 
Bayes) clearance estimates of ibuprofen are regressed independently on each covariatc. 



covariate 


linear vs. constant 


spline vs. linear 


ARSS 


P 


ARSS p 


WT 


-3.89 


<0.01 


-0.10 057 


AGE 


-3.88 


<0.01 


-0.13 0.49 


HT 


-3.64 


<0.01 


-0.21 033 


FED 


-0.61 


003 




DRG 


-058 


0.03 




RACE 


-0.46 


0.06 




LOC 


-0.07 


0.45 




SEX 


-0.03 


0.60 





INITIAL SCREENING 

As a preliminary screening step before the GAM analysis, the individual PK 
parameter estimates are regressed independently on each covariate. This step is 
equivalent to the method proposed by Mairre et al. [8] and may give a first impression 
of the relative v importance of each covariate and of the shape of the relationship 
between it and the model parameters. The results of this screening step for ibuprofen, 
with respect only to the parameter clearance, are summarized in Table I. The covariate 
relationships are modeled linearly or with splines. The significance (p) of the linear 
vs constant model and the spline vs linear model is tested using the F test The 
initial screening shows that AGE, WT, HT, DRG, and FED may significantly (p < 
0.05) influence the clearance of ibuprofen and that the relationship between the 
continuous covariates and clearance is best described by a linear model. 



GENERALIZED ADDITIVE MODEL 

The generalized additive model is derived using a step-wise addition/ deletion 
method [12]. This method steps through a series of models according to the following 
algorithm. At each step of the model building each of the covariates is allowed to be 
deleted from the model, or to enter the model in any of several prespecified functional 
representations; e.g., linearly or as a spline. At each step the model is advanced by 
addition, deletion, or replacement of the single additive term that results in the largest 
decrease in the Akaike Information Criterion (AIC). (The AIC is proportional to the 
residual sum of squares from the GAM fit plus a penalty, proportional to the number 
of parameters in the model.) The search stops when the AIC has reached a minimum 
value. We estimate the GAM using the statistical programming environment S-plus 
(version 3.0, Statistical Sciences Inc., Seattle, WA). However, a FORTRAN program 
to perform the GAM analysis is also available (ref. 11, page 307). 

A GAM was derived for all PK parameters of ibuprofen; CI, V, and ka, however 
only the results for CI are shown as an example. For ibuprofen the GAM analysis 
indicated that CI is a linear function of WT and a function of the categorical covariates 
DRG and RACE. The GAM was derived using both a forward search, starting with 
the NULL model and a backward search, starting with the FULL model, which 
includes all covariates. A total of 71 different models were tested for CI. The path 
taken to the final GAM is summarized in Table II. This table also shows some of the 
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Table II. The top part of the table shows the path taken to the final GAM for the clearance 
(CI) of ibuprofen, starting with the NULL model. In each step the term that results in 
the largest decrease in the AIC is added (+) to the model. The bottom part of this table 
shows some models close to the final GAM (model with smallest AIC value). The - in 
the notation of the models means "is a function of" and the + signs are to be interpreted 
figuratively and not literally. 



step 


term 


ARSS 


RSS 


AIC 
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constant 


12.16 


2343 
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+YVT 


-3.89 


8.26 


2005 
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-036 


7.90 


198.2 
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1963 



Summary of model close to the minimum AIC model 



Equation 


AIC 


Cl-WT+DRC+RACE 


1963 


CI-WT+DRC+RACE+LOC 


1973 


Cl-WT+DRC+RACE+SEX 


1973 



models that (based on the AIC) are very close to the model that best describes the 
data. Table II shows that WT has the greatest influence on CI, explaining 32% of the 
observed variability in CI. DRG and RACE are not as important, and only result in 
a small additional change of the residual sum of squares, explaining an additional 
6% of the variability in CI. HT and AGE which appeared to be important covariates 
according to the initial screening do not appear in the GAM. The reason for this is 
that these two covariates are highly correlated with WT (r > 0.940). 

The relationships between covariates and ibuprofen CI for the final GAM are 
shown in Fig. 2. The solid lines displayed in this figure are the additive contributions 
of each of the terms in the GAM, i.e., the </(. ) for each of the covariates. The actual data 
plotted represents the covariate vs partial residual values. The latter is the individual 
CI estimates minus the predictions based on all other covariates apart from the one 
shown. 



NONMEM ANALYSIS 

In the final (NONMEM) step of the model-building, the non-linear mixed effect 
population model is derived, based on the regression models found in the GAM step. 
Two strategies (roughly corresponding to backward elimination and forward addi- 
tion) can be used in the NONMEM step: (i) include in the initial NONMEM model 
all covariates (with associated functional relationships) that appear in the GAM with 
the lowest AIC value, and in models close to it, and attempt, thereafter, to eliminate 
redundant covariates; or (ii) start with the final GAM as the initial NONMEM model 
and attempt, thereafter, to add non-redundant covariates, selecting them from co- 
variates appearing in models close to the final GAM (see e.g. Table II). The potential 
covariates can be tested in the order of their importance in the GAM analysis. We 
prefer the second approach, especially when the number of candidate covariates is 
large. The parameter estimates of the GAM analysis are used as initial estimates 
for the NONMEM analysis. Model selection is done on the basis of the NONMEM 
objective function, -2 log likelihood (-2LL). The difference in -2LL between a full and 
reduced model is asymptotically \ 2 distributed with degrees of freedom equal to the 
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Fig. 2, Partial residual plot of each of the additive terms of the final GAM for the O of ibuprofen. 
The solid line in each plot shows the contribution of a term to the additive predictor of CI. The actual 
data in each of the plots are the partial residuals; i.e., the observed individual O values minus the 
prediction based on all terms apart from the one shown. Each curve has been centered to have an 
average value of zero. The same scale is used for the y-axis in each plot so that the relative importance 
of the covariates can be compared. The left panel shows the relationship between WT and CI. The 
solid line represents the contribution to the GAM, characterized by a natural cubic spline with two 
internal knots at the 33% and 66% quantiles: 11 kg and 17 kg. The middle panel shows the effect of 
the two dose level groups (DRG, 1= 5 mg/kg and 2= lOmg/kg). The right panel shows the effect of 
race (RACE, l=caucasian, 2=black). For these categorical covariates the points are randomly shifted 
by small amounts on the abscissa to display them better. 



Table m. Change in -2LL when each of the parameters of the covariates appearing in 
the final NONMEM model for ibuprofen is set to zero. 

Covariate WT DRG RACE 

-2ALL(vs.fl=0) 115.0 16.0 9.9 



difference in number of parameters between the two models. Covariates are added 
to the model if they significantly decrease the -2LL ( > 6.6; p < .01), or deleted from 
the model if -2LL does not increase significantly. Covariates are also deleted from 
the model if ± 1 (asymptotic) SE of the parameter estimate includes zero or when the 
magnitude of the covariate effect is small on clinical grounds. 

The NONMEM step in the ibuprofen model building exercise results in the 
following population model for CI 

CI = (0i + 02 (WT - 15) + 0 3 DRG + QaRACE) • expfof) (7) 

Table III summarizes the change in -2LL if each of the 0, in the final NONMEM 
model is set to zero, showing the significance of each of the parameters. Including 
the final set of covariates in the NONMEM model resulted in a decrease of the inter- 
individual variability from 46% to 29% for CI. No additional covariate results in a 
significant change in -2LL. 

Table III also shows that the GAM analysis not only provided all the important 
covariates for the final NONMEM model, but also that it correctly identified the 
relative (order of) importance of each (compare Table III and Table II), and the correct 
shape of the relationship. To check the latter, the relationship of WT to ibuprofen 
CI (the most important relationship) was modeled in NONMEM with a spline. This 
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Table IV. Summary of covariates found to significantly influence clearance according 
to the GAM and NONMEM analysis for all four drugs. The covariates are listed in their 
order of appearance in the GAM. The total number of covariates evaluated was 8 for 
ibuprofen, 9 for quinidine, 13 for prazaosin, and 16 for pcfloxacin. 
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nonparametric representation did not result in a significant decrease in -2LL, con- 
firming the results from the GAM analysis that this relationship is indistinguishable 
from a linear one. 

Similar results were found for V. The GAM step provided three covariates (WT, 
LOC, and DRG) and their functional relationships, all of which are retained in the 
final NONMEM model. One additional covariate, SEX, that is not present in the 
final GAM model was found to significantly influence the fit in the NONMEM step. 
However, this covariate does appear in the GAM model closest to the minimum AIC 
model. 



VALUE OF THE GAM STEP 

To decide if the GAM step is worth the trouble, one must consider whether it 
does what it is designed to do, i.e., identify the important covariates and the shape 
of their functional relationships to parameters; and whether it does so better than a 
simpler approach. Our impression is that the answer to both questions is "yes". 

Regarding achieving its goal, Table IV compares the covariates that are selected 
by the GAM analysis to significantly affect the CI of the four drugs and the covariates 
that are significant in the NONMEM analysis. For all four, the GAM step either 
selects all the important covariates or includes them in GAMs close to the minimum 
AIC model. The latter case happened only once. For quinidine, RF appeared in a 
model close to the minimum model. The GAM also identifies the relative importance 
of the covariates. The covariates that provided the greatest decrease in AIC during 
the GAM analysis also resulted in the largest change in -2LL in the NONMEM 
step. Conversely, those covariates that are of borderline significance in the GAM 
analysis are of minor significance in the final NONMEM model. Similar results 
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Fig. 3. Partial residual plots for the contribution of HT and AGE to the final GAM for the clearance 
of prazosin (upper panels) and for the contribution of HT and AGE to the additive NONMEM model 
(lower panels). The relationships are described by a natural cubic spline with two internal knots at 
the 33% and 66% quantiles: 168 cm and 177 cm for HT and 52 and 58 years for age. Note that the 
GAM plots are shifted on the y-axis in comparison to the NONMEM plots because they are centered 
to have an average of zero. The NONMEM plots show zero effect at the lowest values for HT and 
AGE. See legend to Figure 2 for further details. 



were found for V. So far we have observed that the GAM analysis is inclusive. 
For all four drugs a total of 27 covariates were found to be of significance in the 
NONMEM models. Of these covariates, 25 were found in the final GAM. The other 
2 were found in models very close to the minimum. Regarding the shape of the 
functional relationships, prazosin and pefloxacin provide examples. For the former, 
the GAM analysis identifies important nonlinearities in the relationship between CI 
and the covariates HT and AGE. The spline contribution of each of these covariates 
to the GAM is shown in Fig. 3. A similar nonlinear relationship was included 
in the NONMEM model, yielding a significant decrease (30.0 points) in -2LL (4 
additional parameters) when compared to a linear model for both covariates. The 
spline relationships found in the NONMEM step are also shown in Fig. 3. Note 
the similarity between the splines found in the GAM and NONMEM steps. For 
pefloxacin the GAM analysis indicates a nonlinear relationship between CI and WT 
A similar relationship is found in the NONMEM step, resulting in a significantly 
better fit than with the linear model. 

Regarding whether the GAM analysis provides a better starting point than a 
simpler approach, we compared it to the one-covariate-at-a-time initial screening 
step. For pefloxacin 13 out of the 16 covariates are found to significantly affect CI 
according to the one-covariate-at-a-time analysis (p < 0.05). The GAM analysis 
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Fig. 4. Plot of empirical Bayes estimates of the individual clearance (CI) values for ibuprofen obtained 
from the basic PK-model without any covariates (x-axis) and the final population model including 
all covariate effects (y-axis). The solid line is the line of unity. 

narrows this down to 8 covariates of which 6 are finally selected by the NONMEM 
analysis. One of the covariates that is selected by the NONMEM analysis, SPB, is 
not selected by the one-covariate-at-a-time analysis, but does appear in the GAM. 
Similar results are found for the other drugs. In general, the one-covariate-at-a- 
time analysis tends to select numerous covariates that are not significant in the final 
NONMEM analysis, probably because they are collinear with one or a combination 
of other covariates. Of greater importance, perhaps, is that the one-covariate-at-a- 
time analysis disregards some of the covariates that appear to be of importance in 
the final NONMEM model. The reason may be that their effect is masked by some 
other covariate or by collinearities among covariates. 

These results show that the GAM analysis is a robust approach to select a good 
initial model for the NONMEM analysis. The stepwise addition/deletion method to 
derive the GAM model is considerably quicker than a similar model search within 
NONMEM and provides a nice graphical representation of the relationships between 
covariates and model parameters. 



DISCREPANCIES BETWEEN GAM AND NONMEM STEPS 

There are several factors that may influence the success of the GAM step: (i) It 
may depend on the quality of the empirical Bayes estimates of the PK-PD parame- 
ters. These estimates tend to be "shrunken" towards the population mean, and may 
therefore, fail to reveal important relationships. To the extent that this problem can be 
detected in the NONMEM step, it appears not to be present in the data sets analyzed 
here. One would predict that it might appear when individual data are particularly 
imprecise or few in number, (ii) The empirical Bayes estimates will change as the 
model is built, because the prior changes. For the four examples presented here, we 
have not found this to be a problem. Figure 4 shows that the individual posterior 
Bayes estimates of ibuprofen CI change only marginally when the prior distribution 
is changed from that of the basic NONMEM fit (just the structural model without 
covariates) to that of the final NONMEM fit (including all covariates). As for the pre- 
vious problem, this might prove more important when individual data are imprecise 
or sparse, (iii) The individual PK-PD parameter estimates are weighted equally in 
the GAM analysis, but not in the NONMEM analysis. This could be solved by using 
the precision of the empirical Bayes estimates as weights in the GAM analysis. How- 
ever, we have not found this to be a problem; (iv) The covariates may change within 
one individual from day to day, thereby obscuring covariate-parameter relationships 
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bv "averaging" them within individuals. This can be solved by dividing each indi- 
vidual's data into disjoint contiguous time periods in each of which the covariates 
are fairly constant. We have applied this technique successfully to the prazosin and 
pefloxacin data. Clearly, however, the diff.culty of trying to detect within-.ndiv.dual 
correlation between parameters and covariates will be exacerbated as before if data 
are imprecise or sparse, as now one requires sufficient data in each period to allow 
a reasonably stable empirical Bayes estimate; (v) The AIC is used as ^selection cri- 
terion in the GAM analysis, whereas the log likelihood test is used in the NONMEM 
analysis. The AIC tends to select larger models under the operating conditions o f our 
NONMEM step; hence, relative to the latter, the GAM analysis is inclusive, which is 
an advantage for a screening strategy. 

INTERACTION BETWEEN COVARIATES 

A disadvantage of the additive structure of the GAM is that interactions between 
covariates are not included. To take interactions into account, the ge ™ral s^cture 
of a multidimensional smooth function as described by equation (4) could be used. 
However, this will result in the identifiability problems discussed earher Recently, 
a new approach called PI has become available for covariate selection [14]. This 
approach automatically selects combinations of covariates showing important inter- 
actions and determines their functional relationship using splines. The PI method 
models the interaction between covariates according to the following equation. 



M N 

<t>i = E II 2m.n <■*«> + "» 

m=l n=l 



(8) 



where g m „(.) is a smooth univariate function, in this case a spline. We have used this 
tecLq^; which is implemented in the FORTRAN program PIMPI^ [14], I »sa*en 
for interactions between covariates in the four data sets, as follows: (i) In a first step 
empirical Bayes estimates of m were obtained using as prior the model obtained in the 
NONMEM step. These* represent the residual variability in the PK-PD parameters 
unexplained by the current additive model; (ii) Thebestfit of all possible P^*£f 
two covariates (N in equation 8 equals 2) to the Bayes estimate of * foreach PK-PD 
Parameter is obtained'using the PI method; (iii) The previous NONMEM model is 
updated by including the candidate interaction determined in the step (n). The pair 
of steps (iiHiii) can be repeated, first for other two-way interactions and then for 
three-way interactions, etc., until no interaction further improves the fit. 

For all four drugs, no interactions of this type were found. This suggest; that an 
additive model is often a useful approximation to the more complex true multivariate 
regression surface. A possible reason is that our data are often not good enough (too 
noisy and too sparse) to require complex interactions for adequate description. 

The PI method can only be used to screen for interactions between continuous 
covariates. Interactions between categorical covariates and categorical and contin- 
uous covariates are therefore analyzed using S-plus. For prazosin an interesting 
interaction was found. The final additive model in NONMEM showed that Q* a 
function of HT, AGE, HCTZ, and RACE and that V is a function of HT and HCTZ. 
An analysis of the empirical Bayes estimates of * using this additive model as a pnor 
showed a profound interaction between HCTZ and SEX for both CI and V Figure 5 
shows the interaction between SEX and HCTZ for the CI of prazosin. Add-on of 
a saturated interaction model between HCTZ and SEX to the (additive) NONMEM 



O 
TJ 

(0 
0> 



o 

CL 



20 n 
15 

10 A 
5 
0 
-5 
-10 -| 
-15 



Fig. 5. Partial residual plot showin 
prazosin. The partial residuals are pi 
which indicates that the individuals 
mean partial residual values for each. 

model results in a significant d 
This model can be further refir 
not significant and by deleting 
resulting, finally, in a model wi 
with one less parameter. Thu 
NONMEM model results in th 
that appear to be significant in 
probably substituting for the ir 

DIAGNOSTICS FOR NONLI 

The estimation methods ir 
sions of least-squares method: 
sensitive to "outliers". Thus, r 
only a few influential observati 
hires of those observations rati 
and parameters. It is therefore \ 
an unusual large influence on t 
focuses on the relationships bet 
cific covariates, influence is na 
than the level of single observa 

In regression analysis, cas< 
detect influential observations 
propose to use two diagnostics 
gression analysis: the Cook-sco 
case-deletion diagnostics, as w 
the data of one individual on i 
overall influence of the deletion 
0 now is used generically to re 
and random effect, of the popv 
equations: 

COOKi = (\ 



80 



solved by dividing each indi- 
i each of which the covariates 
iccessfully to the prazosin and 
ng to detect within-individual 
e exacerbated as before if data 

data in each "period" to allow 
AIC is used as the selection cri- 
i test is used in the NONMEM 
the operating conditions of our 
1 analysis is inclusive, which is 



AM is that interactions between 
> account, the general structure 
by equation (4) could be used, 
ms discussed earlier. Recently, 

covariate selection [14]. This 
riates showing important inter- 
• using splines. The PI method 

to the following equation. 



case a spline. We have used this 
urogram PIMPLE [14], to screen 
ets, as follows: (i) In a first step, 
; prior the model obtained in the 
ability in the PK-PD parameters, 
»est fit of all possible products of 
?s estimate of 77, for each PK-PD 
2 previous NONMEM model is 
rmined in the step (ii). The pair 
-way interactions, and then for 
her improves the fit. 
ere found. This suggests that an 
more complex true multivariate 
1 are often not good enough (too 
s for adequate description, 
nteractions between continuous 
ites and categorical and contin- 
js. For prazosin an interesting 
NONMEM showed that CI is a 
is a function of HT and HCTZ. 
ng this additive model as a prior 
SEX for both CI and V. Figure 5 
the CI of prazosin. Addition of 
;EX to the (additive) NONMEM 



1 



o 
en 



O 
Q_ 



20 -| 
15 - 
10 - 

5 - 

0 
-5 

-10 H 



-15 



0 (P 
0 



0 

Q 



10 



1 



male 



Sex 



female 



Fig. 5. Partial residual plot showing the interaction between SEX and HCTZ for the clearance of 
prazosin. The partial residuals arc plotted as "0" which indicates no co-treatment with HCTZ or "1" 
which indicates that the individuals received co-treatment with HCTZ. The solid lines represent the 
mean partial residual values for each of the four groups. 



model results in a significant drop in -2LL of 21.4 on four extra degrees of freedom. 
This model can be further refined by deleting some of the interaction terms that are 
not significant and by deleting RACE from the CL model and HT from the V model, 
resulting, finally, in a model with -2LL 14 points lower than the additive model and 
with one less parameter. Thus, inclusion of the important interaction term in the 
NONMEM model results in the removal of two factors, RACE for CI and HT for V, 
that appear to be significant in the additive NONMEM model. These two factors are 
probably substituting for the interaction effect. 



DIAGNOSTICS FOR NONLINEAR MIXED EFFECTS MODELS 

The estimation methods in NONMEM (and GAM) can be thought of as exten- 
sions of least-squares methods. It is well known that such methods can be very 
sensitive to "outliers". Thus, model building results can be strongly influenced by 
only a few influential observations, and the fitted model may reflect the unusual fea- 
tures of those observations rather than the general relationship between covariates 
and parameters. It is therefore very important to detect those observations that have 
an unusual large influence on the results of the analysis. Since population analysis 
focuses on the relationships between individual PK-PD parameters and patient spe- 
cific covariates, influence is naturally defined at the level of the individual, rather 
than the level of single observations. 

In regression analysis, case-deletion diagnostics are the standard tool used to 
detect influential observations [15,16]. For a non-linear mixed effects model we 
propose to use two diagnostics that are generalizations of those used in standard re- 
gression analysis: the Cook-score (COOK) and covariance-ratio (COV-ratio). These 
case-deletion diagnostics, as we use them, measure the influence of the removal of 
the data of one individual on the parameter estimates of the current model. The 
overall influence of the deletion of the ith individual on the estimate of 8 (the symbol 
0 now is used generically to refer to the vector of all parameters, both fixed effect 
and random effect, of the population model) can be calculated from the following 
equations: 



COOIU = {y/(0i - OfcoviOyHei - 6)) (9) 
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det(c<w(0,)) 



det (cov(0)) 



(10) 



where 0, is the vector of parameter estimates when the tth individual is deleted, 
cov{.) denotes the covariance matrix of the parameter estimates, and det{.) denotes 
the determinant of a matrix. The Cook-score for the ith individual measures (approx- 
imately) the average absolute change in the value of the parameter estimates when 
the ith individual's data are deleted, each parameter change scaled to its uncertainty 
(standard error) when estimated using all the data. The COV-ratio for the ith indi- 
vidual measures the "overall" uncertainty (standard error) of all parameter estimates 
on deletion of the ith individual's data, scaled to the full-data overall uncertainty. 
A large Cook-score indicates that the individual's data have an inordinate effect on 
parameter estimate values, and hence are to some degree inconsistent with the rest 
of the data, while a large value of the COV-ratio indicates that the individual's data 
contribute considerable information to the fit (are influential). A small value of the 
Cook-score indicates that the individual's data are compatible with the rest of the 
data, while a small COV-ratio indicates that parameter precision improves when the 
individual's data are deleted, suggesting that they are incompatible with the other 
data. 

To detect the influence on a single parameter 0 k of the model, similar diagnostics 
can be calculated using the following equations 



COO!u kJ = 



(ft,. - Ok) 2 

var (Ok) 



COV-ratio^.. = 



var (0 k j) 



\ var(6 k ) 



01) 



(12) 



where var(.) denotes the variance of the parameter estimate, and 0*,. denotes the 
estimate of 0 k when the data of the ith individual are deleted. The Cook-score 
for a single parameter and tth individual measures the change in the value of the 
parameter estimate when the ith individual is deleted, scaled to the uncertainty 
(standard error) when the parameter is estimated using all the data. The COV- 
ratio for a single parameter and ith individual measures the uncertainty (standard 
error) of the parameter estimate on deletion of the ith individual's data, scaled to the 
uncertainty when estimated from the full-data. The implications of large and small 
values of these individual parameter diagnostics for individual subjects' data are the 
same as for the overall diagnostics, except that they are specific for the particular 
parameter. . . 

The diagnostics presented here are not intended to provide rules for the rejection 
of data. Influential subjects' data should not necessarily be eliminated from the 
analysis. These individuals' data may simply be particularly informative. However, 
influential subjects may point out weaknesses in the analysis and may suggest that 
the current model is inadequate or may suggest that additional data need to be 
collected. For example, if the removal of an influential individual changes the value 
of an estimated parameter to such an extent that its effect becomes insignificant, 
one is uncomfortable extrapolating the full-data parameter value to future subjects 
as its value depends, essentially, on just one individual. In general, influential 
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Fig. 6. Plot of the overall Cook-score and COV-ratio when each of the individuals is deleted from 
the NONMEM model for ibuprofen. The numbers plotted are the actual index number for each of the 
individuals (up to 92). 

individuals are those whose data are far removed from other individuals in the 
(high dimensional) observation-space; identification of regions of observation space 
with inadequate coverage for reliable estimation and prediction is an important data 
analysis goal [15]. 

In standard regression analysis, numerical methods are available that calculate 
case-deletion diagnostics without the need to re-estimate parameters after removal 
of each observation (individual). For nonlinear mixed effect models, however, no 
such methods are currently validated and case-deletion diagnostics are best calcu- 
lated by brute force, i.e., the model is refit with each individual deleted. We have 
developed a NONMEM control stream that automatically refits the model with each 
of the individuals deleted and that saves the resulting case-deletion parameter es- 
timates and corresponding variance-covariance matrix of the estimates [17]. The 
case-deletion diagnostics can then be calculated using Eqs. (9M12). 

EXAMPLE OF CASE-DELETION DIAGNOSTICS 

We here present some results of using the case-deletion diagnostics described 
above to examine the ibuprofen analysis, and also some results of so doing for the 
other 3 drugs. 

According to the NONMEM step of the ibuprofen analysis, CI is a function of 
WT, DRG and RACE (Eq. 7), and V is a function of WT, LOC, DRG, and SEX: 

V = (0 5 -f MWT-15) + (hLOC + e 8 DRG + 0 9 SEX)-exp(r) v ) (13) 

This model was refit to the data with each of the individuals deleted. The values of 
the overall Cook-score and COV-ratio are shown in Fig. 6 for each individual. Figure 
6 identifies certain suspect individuals (those appearing along its "edges"): number 
84 has a particularly low COV-ratio (outlier?), number 22 has a high COV-ratio 
(influential?), and numbers 31 and 47 have high Cook-scores (outlier? influential? 
both?). These individuals can be further evaluated by plotting the Cook-score vs 
COV-ratio for each of the individual parameter estimates. Figure 7 shows these 
plots for the effect of WT (0 6 ), LOC (0?), DRG (0 8 ), and SEX (0 9 ) on V. These figures 
clearly indicate that the data of individuals 31, 84, and 47 have an unusually great 
effect on the value of the estimate of V when compared to other individuals. For 
example, individual 47's data changes 0 9 almost 1.5 times the standard error of the 
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Fig. 7. Plot of the Cook-score and COV-ratio for the parameters describing the effect of WT <0 6 ), LOC 
(b), DRG ($ s ), and SEX (0g) on the volume of distribution of ibuprofen. The numbers plotted are the 
actual index number for each of the individuals (up to 92). 

parameter estimate. Table V summarizes the parameter estimates for the model 
for V when each of these individuals is deleted from the analysis. Removal of 
individual 31's data renders the effect of DRG insignificant, whereas removal of 
individual 47's data renders the effect of SEX insignificant Both effects are highly 
significant (p<0.01) according to the likelihood ratio test when all individuals are 
included in the analysis. Furthermore, the magnitude of these (binary) effects in 
the full-data analysis is notable: both change the average V by about 18%. The 
number of males/females was 53/40, while the number of subjects with low/high 
DRG was 45/48. Since all these numbers are appreciable, the fact that the apparent 
effects of the covariates change so profoundly upon deletion of just one individual 
in each case suggests that the "covariate effect" is not real, but simply an artefact 
due to outlying individuals. If, however, there had been severe covariate imbalance; 
i.e., only a few individuals with the same SEX value as number 47, or only a few 
individuals with the same DRG value as number 31, then one might suspect that the 
effect could be real, and seek confirmation through study of additional individuals 
with like covariate values. The effect of removal of the data of individual 84 on 67 
is also profound. LOC is a highly significant factor when the data of all individuals 
are analyzed, changing -2LL by 50 points when included in the model. Yet one 
must conclude that the magnitude of its effect is poorly defined, since it changes by 
almost 1 standard error when subject 84 is removed. In all but the last of the above 
cases, it is true that the standard errors of the parameter estimates based on the full 
data are sufficiently large that the covariate effects are suspect, but the case-deletion 
diagnostics add considerable specificity to this vague suspicion. 

Case-deletion diagnostics will identify influential individuals even though they 
be otherwise difficult to identify due to the complexity of the population data. Even 
may plots of the empirical Bayes estimates of the model parameters vs the covariates 
such as-in Fig. 2 may not clearly identify influential subjects because of the complex 
influence structure of each subject on all parameter estimates. 

Case-deletion diagnostics are also helpful in explaining some of the prazosin 
results. For this drug, several covariates are found to significantly influence CI (see 
Table IV) and V, resulting in a marked drop in -2LL when they are included. However, 
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Table V. Effect of removal of certain individuals' data (#31, #47, #84) on the parameter 
estimates of the NONMEM population model for the volume of distribution of ibupro- 
fen compared to the parameter estimates based on all individuals' data (all #). The 
figures in each cell are the parameter estimate (standard error of estimate). The change 
in -2LL (based on the analysis with all #) when each of the parameters estimates (0) is 
set to zero is also shown (first column). 
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the inter-individual variability in CI and V is not reduced by inclusion of these 
explanatory covariates, a perplexing contradiction (unexplained inter-individual 
variability is being "explained", why then does its magnitude not decrease?). Case- 
deletion diagnostics reveal that this is due to one subject. Removal of this subject 
results in a 25% decrease of inter-individual variance in CI and a 50% decrease in the 
standard error of this estimate. Similar results are found for the estimate of inter- 
individual variability in V. This subject did not have unusual CI and V estimates when 
compared to other subjects. However, in contrast to the others, this individual's CI 
and V were not well explained by the final population model. 

CONCLUSIONS 

In this paper we have shown that a preliminary GAM analysis is a robust 
approach to selecting a good initial model for a subsequent NONMEM analysis. 
The GAM step, selects the most important covariates and indicates important non- 
linearities in the relationships between covariates and PK-PD model parameters. The 
method is efficient and provides illustrative graphics. A disadvantage of the additive 
structure of the GAM is that interactions between covariates are not included, but 
we here present a screening method to detect such interactions. When it is applied to 
our four data sets, no interactions are found, suggesting that a GAM may often be an 
adequate approximation to an inevitably more complex true multivariate regression 
surface. 

We have also shown that case-deletion diagnostics can detect individuals whose 
data is unusually influential in the final NONMEM analysis. Use of such diagnostics 
should allow population modeling to become "more robust. 
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